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Using amplitude equations, we show that groups of identical nano-mechanical resonators, inter- 
acting with a common mode of a cavity microwave field, synchronize to form a single mechanical 
mode which couples to the cavity with a strength dependent on the square sum of the individual 
mechanical-microwave couplings. Classically this system is dominated by periodic behaviour which, 
when analyzed using amplitude equations, can be shown to exhibit multi-stability. In contrast 
groups of sufficiently dissimilar nano-mechanical oscillators may lose synchronization and oscillate 
out of phase at significantly higher amplitudes. Further the method by which synchronization is 
lost resembles that for large amplitude forcing which is not of the Kuramoto form. 
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I. INTRODUCTION 

, Synchronisation of coupled oscillators arise in many different contexts in biology, chemistry and engineering [l|, Q • 
Such systems show surprising emergent behaviour and can be used to encode and process information Q. In this 
paper we show how synchronisation can arise in arrays of nano-mechanical resonators interacting via a common 
Q electromagnetic field mode. Recent progress in opto-mechanical and nano-mechanical systems now enables very high 
' frequency mechanical resonators to be coupled strongly to one or more modes of the electromagnetic field in a resonant 
cavity This is largely driven by a desire to explore the deep quantum domain in which the mechanical resonator 

is prepared at or near its vibrational ground state As the coupling is essentially nonlinear, the resulting classical 

dynamics can be complex and must be thoroughly understood if one is to make sense of the quantum phenomenon. 

The common feature in these systems is the so-called 'radiation pressure coupling', whereby the displacement of 
each mechanical resonator independently changes the resonance frequency of a common electromagnetic cavity field 
by an amount proportional to the displacement of each mechanical resonator. This means that there is an effective 
'— — conservative force acting on each mechanical resonator proportional to the circulating power in the electromagnetic 
. cavity. If the cavity is externally driven, this interaction mediates an indirect all-to-all coupling between each of the 



mechanical resonators that is highly nonlinear. 

If the oscillators are identical, a collective variable can be used to understand the dynamics. In this paper each of 
the oscillators is a bulk flexural vibrational mode of a mechanical resonator. The resulting set of equations is similar 
to that considered by Marquardt et al Q , who found that multi-stability was an important feature of the dynamics 
for small mechanical damping. Here we are able to derive amplitude equations for the collective variables and use 
' these to map out regions of multi-stable behaviour in the system. 

For nonidentical phase oscillators Kuramoto [9( used a collective variable (Kuramoto's order parameter) to char- 
acterise the synchronisation between the oscillators. More recently the collective dynamics of opto-mechanical arrays 
has been described by Heinrich et al. [Io| . who give some results on synchronisation based on a phase model related 
to Kuramoto's model. Like our model, this paper is based on the radiation pressure coupling between the field and 
mechanical elements. Unlike our model, the mechanical resonators in Heinrich et al. interact with a local electro- 
magnetic field mode and are directly coupled by elastic forces. The more complex coupling in our model results in 
a different mechanism for the loss of synchronisation which typically occurs for large amplitude forcing not small 
amplitude forcing as occurs in the model of Heinrich et al Nevertheless we are able to give specific results on 
synchronisation for two and three mechanical resonators interacting via a common cavity mode and relate these to 
the behaviour of a collective variable, which is related to the cavity field amplitude. 

Much of the previous work on synchronised nonlinear oscillators is based on a direct, usually nearest-neighbour, 
interaction between the individual oscillators and amplitude equations have been used successfully to analyse the 
dynamics of such models [ill ] . We show that amplitude equation methods can also be applied to understand the 
dynamics of the more complex all-to-all coupling that occurs in our model. Synchronisation in coupled micro- 
electromechanical systems (MEMS) has been described [l2j and observed [l3l |. 

There are at least four kinds of physical implementations of the system discussed here. Firstly, in circuit QED, a 
coplanar microwave cavity contains the electric field which forms the common field mode. Nano-mechanical resonators 
can then be placed so as to form one plate of a capacitor with the central conductor of the microwave cavity thereby 
modulating the microwave cavity frequency [TJ- Secondly, at optical rather than microwave frequencies, an opto- 
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mechanical system can be formed by placing micro- mechanical dielectric membranes inside the optical cavity (15| . 
Thirdly, a toroidal optical whispering gallery mode (WGM) cavity is manufactured on a tapered platform raised off a 
substrate [l6j]. The mechanical vibrations of the toroid modulates the frequency of the WGM. Finally, opto- mechanical 
phononic crystals can be fabricated which are simultaneously photonic crystal lattices to produce localised optical 
and mechanical modes [171, [l8| . 

In the bulk of this paper we will consider a nano-mechanical system where a single mode of a superconducting 
microwave resonator is coupled to the displacements of TV nano-mechanical resonators. However, the dynamical 
model we derive in this section is applicable to the other physical implementations in different experimental contexts. 
In general our model applies to a system of TV + 1 oscillators: N single flexural modes of independent mechanical 
resonators whose displacements are coupled to a common single electromagnetic field mode, also modelled as a single 
simple harmonic oscillator. The coupling between each mechanical resonator and the microwave field in the cavity is 
capacitive and results in a frequency shift of the cavity resonance frequency that, to lowest order, is proportional to 
the displacement of the mechanical resonator. This results in a force on each mechanical resonator that is proportional 
to the intensity of the microwave field in the cavity. This is often called radiation pressure coupling [19]. A schematic 
of this system is given in Fig. [1] 



Ground plate 




FIG. 1: A schematic of the nano-electromechanical system under consideration. A superconducting microwave cavity of 
frequency uj c mediates a coupling between N nano-mechanical resonators capacitively coupled to it. The ith nano-mechanical 
resonator has resonant frequency Ui and microwave-mechanical coupling strength Qi. The microwave cavity is driven by a linear 
drive of amplitude e at a detuning from the cavity of <5. 



We model the dynamics of the microwave field in the co planar transmission line by a lumped circuit LC electrical 
resonator and the dynamics of each mechanical resonator is modelled as a single simple harmonic oscillator. The 
Hamiltonian for a single nanomechanical resonator interacting with the microwave field is 

^ 2 Q 2 / \ s~\ P 2 muj2 ■? ,«S 

U =2L + ik) +V{t)Q + h + — q ' (1) 

where the first term is the inductive energy with the <P the flux through the equivalent inductor with inductance 
L. The second term is the charging energy with Q the charge on the equivalent capacitor with capacitance C(q), 
which varies with the displacement of the mechanical element. The third term represents the potential energy due 
to an external AC bias voltage of the equivalent circuit resonator. The fourth term is the kinetic energy of the 
mechanical resonator of effective mass m and the last term is the elastic potential energy of a single flexural mode 
of the mechanical resonator with u). As the displacement is small compared to the equilibrium distance between the 
mechanical resonator and the central conductor of the microwave cavity we can expand C(q) to linear order in q 
around the equilibrium displacement qo to get an effective Hamiltonian 

n = ^- + $ r + ^+ mu>Y + AQ 2 q + v(t)Q , (2) 
2L ZCq Ira 
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where C = C(qo) and A = -\ 



. The classical Hamilton's equations are 

9=90 

d<P O 

- = % + 2A Qq + v{ t), 
dQ = _£ 

dt L ' /q\ 

dq P (3) 



dt m 
dp 
dt 



-mw q — AQ . 



When A = 0, the circuit equations of motion describe simple harmonic oscillation at the frequency 

1 



o 



(4) 



It is convenient at this point to define dimensionless canonical variables. We do this by first fixing two energy scales, 
one for the circuit degrees of freedom E c and one for the mechanical degrees of freedom, E m . The dimensionless 
canonical variables, (x c ,y c ) for the circuit and (x, y) for the mechanics, are then defined by 



Q 



q (5) 



2E m 
muj 2 



y = 



We now anticipate an eventual quantum mechanical treatment and set E c = ftw Cl E m = hw. The appearance of H at 
this stage does not signify anything more than a convenient conversion factor between energy and frequency. We also 
define a complex amplitude for the circuit degrees of freedom as 

a = x c + \y c , (6) 

in terms of which we can write the Hamilton's equations of motion as 



da 

dt~ 

dx 

dt 
dy 

dt 



—\uj c a — ig (a — a*)x + £{t) , 

uy , (7) 



-ujx - gyl , 



where 



V2AC E C 

9 = 



£{t) 



\frnEn 
v(t) 



(8) 



^/2E~L ' 

We now assume that the circuit is harmonically driven and set 

£{t) = £q sinaj£)t (9) 

and define the rotating variable a = ae lulDt (equivalent to going to the interaction picture in the quantum description). 
If we then drop rapidly rotating terms (compared to the time scale of observations) , the equations of motion may be 
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approximated by 

da . _ 

— — = —Ida — igax — ie , 
dt 

£ = (10) 

d V 5,-,2 

where e = — and the detuning 5 = w c — Up. Noting the the time averaged energy of the energy in an LC circuit 
is proportional to |a| 2 we see that the effective coupling between the microwave field and the mechanical resonators 
is described by the effective Hamiltonian g |a| 2 x. This form of coupling is often termed radiation-pressure coupling 
flTi being proportional to the circulating power in the cavity field. 

We include dissipation of both the microwave field mode and the nano-mechanical resonators using the quantum 
mechanical master equation to incorporate fluctuations correctly. We derive this in section HTT1 There we show that in 
the classical description, the systematic effect of damping (i.e. ignoring fluctuations) change the Hamilton equations 
to 

da . 

— = — id a — wax — ie — Ka , 
dt y 

— =uiy-yx, (11) 

dy g , |2 

—— = — lux — — a — -fy , 
dt 2 1 1 ,y ' 

where k, 7 are the energy decay rates for the electrical and mechanical energy respectively and we have dropped the 
bar as from this point on we simply take it for granted that we are working with the rotating variables for the cavity 
field. 

In this paper we are interested in the dynamics of N mechanical resonators interacting with a single mode of the 
microwave field in the circuit. Assuming a coupling of the form ^Zi9i x i l a | 2 j we see that the equations of motion may 



be expressed in terms of collective variables. 



A M 

—— = —ida — ie — ia > NiXi 
dt f-^ 



Ka , 



d ^ v v (12) 
dt 

dYi Gi . .2 

— = -w,!, - — |a| - jiYi , 



where Xi, Yi and Gi are the M collective combinations 



1 x ^ 



AT, Z-^ 



at ^9 3 Vj, (13) 

1 jest 



Ni 

and Si are collections of Ni identical nano-mechanical oscillators with individual classical positions and momenta 
Xj and yj respectively. We note that the other experimental contexts mentioned in this introduction can also be 
described by the same differential equations (|12[) . For example, multiple opto-mechanical membranes in an optical 
cavity with are described by these equations with different resonant frequencies and coupling strengths [151 ] . We give 
a list of the achievable experimental values for various experiments in table U in section [X] 

In the following section, we present a detailed analysis of the steady state structure of the nonlinear semi-classical 
system, including local and global bifurcations. Since the behaviour is dominated by oscillatory motion, amplitude 
equations are derived from which we can obtain specific results about the existence and stability of periodic orbits. 
It is then a simple step to derive coupled amplitude equations for the case where the mechanical oscillators are not 
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identical and we analyse two and three coupled oscillator systems. In section Hill we give a quantum description of the 
many-body system, and calculate the steady state quantum noise spectra as the first stable limit cycle is approached. 
Finally in section section Hvl we summarise our results and suggest new directions for further work. 



II. DYNAMICS OF THE CLASSICAL MODEL 



Although there are regions of the parameter space where stable critical points exist, periodic motion plays a major 
role in the dynamics for the cases of both the identical and the nonidentical resonators. If the mechanical resonators 
are identical, even if their couplings are nonidentical, they will synchronize, in phase, to form a single mechanical 
mode. However the synchronized motion exhibits multi stable behaviour. The first two sections, below, discuss the 
synchronized motion of identical mechanical resonators (|14[) . largely via amplitude equations. If, on the other hand, 
the mechanical resonators naturally oscillate at different frequencies, desynchronization can occur. To analyze this 
we consider the synchronization between different frequency groups. The resonators can then be attracted to out-of- 
phase solutions that oscillate at much greater amplitudes. In the final section we obtain specific results, via coupled 
amplitude equations, for synchronization between two and three frequency groups. 

For all of the bifurcations that occur a scaled version of the cavity forcing e, which is tunable in an experiment, 
can be thought of as the bifurcation parameter. There are two time scales in the system; the amplitude decay rate 
K of the common cavity mode and the decay rate of the resonators, which is an order of magnitude smaller and will 
be important for the derivation of the amplitude equations. The amplitude decay rate k of the common cavity mode 
provides a natural time-scale and we introduce: a new time parameter t' = ftt; re-scaled nano-mechanical variables 
X[ = ^ and Y( = and dimensionless coupling constants 5' = ~, e' = ^, = ^r, 7- = ^f, G- = and 

= \l + 7; 2 . This gives 

a M 
da 



<=i (14) 

If the uncoupled mechanical resonators are identical = u>', j[ — 7' =>■ uli — a)'), then the oscillators synchronize. 
This is a natural consequence of linear damping and the fact that each oscillator experiences the same forcing. Consider 
u = X[ — Xj then u = is a stable solution of its equation of motion: 

d 2 u _, 2 „ ,du 



■u> u 



27'^, (15) 



dt' z dt 
provided 7' > 0. 

The synchronized motion can then be represented in collective variables (|14[) which, suppressing the use of primes, 
gives the following 



da 
~dt 

d 2 X Gu) , ,2 „ dX 



(1 +16) a- laNX - ie , 

(16) 



, = -u z X - — \a\ 2 - 2 7 — . 

dt 2 2 1 1 1 dt 

For the remainder of this paper we suppress the uses of primes in the notation, and remind the reader that all couplings 
are now dimensionless with the cavity deca y rate determining the natural time-scale of the system. 

From a dynamical point of view eV NG acts as one parameter and in fact both N and G could be removed by 
scaling 



(17) 




So if the number of resonators is increased, smaller values of the driving are necessary to achieve the same effect. 
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A. Critical points, bifurcations and stability 



Without forcing, e = 0, the origin is a stable critical point. As the e is increased from zero the critical point moves 
away from the origin, its position given by the single real root of the cubic 

2uj 2 X (l + (8 + iVXo) 2 ) + Guje 2 = , (18) 

where ao = — \^s+nx ) • However it loses stability on a Hopf bifurcation, creating a periodic orbit, for both 5 > 
(red detuning) and 5 < (blue detuning) provided 7 > and small. The dynamics of this periodic motion is the 
subject of the next section. 

For 8 > (red detuning) the Hopf curve is a perturbation of that from the 7 = case where V NG e = \/28w. To 
first order in 7 it is given by 



e = e H (w, 6, %NG) 
For 8 < (blue detuning) e is order ^/j: 



\ 



2w ( (l+a; 2 r 
NG\ +1 28lo 2 



(19) 



e H (w,S,j,NG) 



\ 



7 (1 + S 2 ) ((<5 2 - uj 2 + l) 2 + 4w 2 ) 



(20) 



-SNGuj 



The Hopf bifurcation is subcritical for 8 < — yj 8 "^+ 3 (blue detuning), where periodic orbits can exist for e < 

€h (w, 5,7, NG). In fact many stable limit cycles can exist for some parameter values because of the presence of 
saddle node bifurcations of limit cycles each creating a stable and unstable pair of limit cycles. This leads to multi 
stable behaviour that has been noticed elsewhere [1, HO] for similar systems. These bifurcations are shown in fig. [2]for 
lj = 2 and 7 = 0.001. The limit cycle bifurcations were produced using using the amplitude equations described in 
the next section, however similar results can be produced by following the limit cycles numerically using the package 
MatCont [2l|. For 8 < (blue detuning) eight of the saddle node bifurcations of limit cycles are shown indicating 
regions where there are 1-8 pairs of stable and unstable periodic orbits. See the caption for specific details. MatCont 
indicates the situation is dynamically more complicated for 8 > (red detuning) involving period doubling and regions 
of chaos. 

Although most of this paper is devoted to the case of blue detuning, where 8 < 0, it is worth mentioning that for 
8 > v3 there is a region where three critical points exist given by the roots of the cubic given above. This triangular 
shaped region 



2u> (28 + VS 2 



sVng^\/s + Vs 2 ^ 



< e < e,„- = 



2w (28 - VS 2 ^) 



sVng^^/s - Vs 2 ^ 



(21) 



is bounded by (e = e sn ±) saddle node bifurcations, shown as green lines in fig. [2j These intersect in a cusp bifurcation 
at 8 = V3 and eVNG^J = 

The case ui — 2 is relevant for the experiments described in 0, However for ui > 2, as in 0, [23[, there is 



no qualitative change in the bifurcation diagram, although the generalized Hopf bifurcation (8 = — 



) occurs 



for larger values of \8\. Fig. [3] shows the corresponding situation for a) u — 5 and b) oj = 10 and 7 = 0.001 with 
8 < (blue detuning) . Multistable behaviour due to the presence of limit cycles stacked above each other remains an 
important feature (see also fig. HJ). 



B. Amplitude Equations and multi-stability for blue detuning (5 < 0) 

Periodic orbits and multiple periodic orbits can exist, if the weakly forced oscillators are sufficiently weakly damped. 
This multi stable behaviour, resulting from the play off between weak damping and cavity forcing, has been noted 
elsewhere [1, HI} - Here we explore it in more depth using amplitude equations. 

The method relies on defining a slow time which is proportional to the damping rate of the resonators: (r = jt) and 
on assuming that the forcing is on the order of the square root of the damping: e = ^p^Z. Then the cavity amplitude 
is naturally of the same order as the forcing and we can obtain equations for the slowly varying amplitude A(t). Let 



X = X + (A{T)e iQt + c.c.) = Ao + 2 \A(t)\ cos (uit + 6) , 



(22) 



7 




FIG. 2: The bifurcation diagram for u> = 2 and 7 = 0.001. In the shaded region there are no periodic orbits and there is one 
stable critical point. The Hopf bifurcation curve, which is in red, provides a partial boundary of this region. At the generalized 
Hopf GH (which is at 8 = \fl for co = 2) the Hopf bifurcation changes from super to sub critical. For 8 < s/l the Hopf 
bifurcation is sub critical and periodic orbits exist to the left of the Hopf curve. Also for < 8 < v3 there are regions where 
periodic orbits exist to the left of the Hopf curve. The blue curves A-GH, BGK, CFK , DEK, KHCusp, KMcusp etc are saddle 
node bifurcations of periodic orbits creating a stable and an unstable periodic orbit existing to their right. (Only the first 8 
are shown.) The lozenge like dashed curves are also saddle node bifurcations of periodic orbits, this time destroying a stable 
and an unstable periodic orbit. (Once again only a sample are shown.) In region ABG(GH) and HM Cusp there is one stable 
critical point and a pair of periodic orbits with opposite stability. In region G(GH)HK there is one unstable critical point and 
one stable periodic orbit. In region BCFG and the region to the left of M Cusp there is one stable critical point and two pairs 
of periodic orbits with opposite stability. In region FGK there is one unstable critical point and two stable periodic orbits and 
one unstable periodic orbit. In region CDEF there is one stable critical point and three stable and three pairs of periodic orbits 
with opposite stability etc. 



where Xq is the critical point of the system given in the previous section, which is 0(e 2 ). Given that 7 is small and 
both e and \a\ are 0(y/j) then X + uj 2 (X - X ) 2jiuj^e iQt + c.c. [H]. The cavity forcing (^ |a| 2 ) can then be 
written as a sum of products of Bessel functions. To see this substitute X = Xq + 2 \A\ cos (uit + 8) into the cavity 
equation; 

— = - (1 + i (5 + NX + 2N \A\ cos (uit + 9))) a-ie. (23) 
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a) w=5 b) 00= 10 




FIG. 3: Bifurcation diagrams a) ui — 5 and b) w = 10 and 7 = 0.001 with S < (blue detuning) showing the Hopf bifurcation 
(red) and saddle node bifurcations of periodic orbits. The labeling in 3a) is similar to fig. [2] For instance in the region 
ABG(GH) there is one stable critical point and a pair of periodic orbits with opposite stability. 



Then if 

a = e i ^Y,B m e im ^ t (24) 

m 

it follows that 

a = -4{t)a + e 1 ^ ^ imwB m e imu ' , (25) 

m 

and using the Jacobi Anger expansion Q X)^L-oo 'i n Jn{z)e mB = e lzcosff , this can be matched to the right hand side 
of the cavity equation if 



■m+l T f 2N\A\\ 

• 2N\A\ , 1 + e M - J 

i/)(t) = — cos (wt + 6) and B m = ^ '- . 

uj k + imw 

where k = 1 + i(6 + NXq) and J m (x) are Bessel functions of the first kind. Substituting this back into 



gives an amplitude equation for the oscillation in terms of sums of pairs of Bessel functions, 



Identical mechanical resonators synchronize to oscillate with amplitude A(r) given by this equation. 
In polar form (A = re 10 ) the equations become 

dr _ 2 ^ .- . [2Nr\ [2Nr\ 

"j- = -r + Ge a mr {5,U>) Jm I ) Jm + 1 [ J ' 



m=0 



d9 Ge 2 ^ ,- . T 2Nr\ 2Nr 
-T- = H — — 2^ tt "» {0, w) J m I — — I J rn+ l I — — 

' ' m=0 \ ' \ 



(26) 



X = -u?X-^-\af -2X (27) 



_ T f 2N\A\\ T { 2N\A\\ 

dA = _ A _iGeV^ ~ M^JjWjj^j _ (28) 
dr 4 ^— ^ (k + i (m + 1) a;) (k* — ima;) 

m— — oo v 



(29) 
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where 



a mr (8,oj) 



a mi (8, uj) = 



8ui 2 (2m + 1) (l + S 2 + ui 2 m (m + 1)) 



((l + P - m 2 uj 2 ) 2 + 4to 2 u; 2 ) ^(l + P - (m + if uj 2 ^ + 4 (m + l) 2 uj 2 ^j 
8u> (2m + 1) ((1 + I 2 - to 2 w 2 ) (l + S 2 - (to + l) 2 uj 2 } + 4m (m + 1) ui 2 ^ 
2 ((l + P - m 2 uo 2 ) 2 + 4mV) ^(l + 6 2 - (m + l) 2 u 2 ^ + 4 (to + l) 2 w 5 



(30) 



7 ((<5 2 -c J 2 +k 2 ) 2 +4k 2 W 2 ) 

and 8 = 5 + Xq. For <5 < (blue detuning) then Xq k, — ^ 2«Zi 2 5 • Since each term in the sum has |A| as 

a factor, the amplitude equation may be rewritten as 

A A 

—— = —A + Ge 2 NAF(N \A\ , u, 8) , (31) 
dr 

where F (Nr,ui, 8) is a complex function. The conditions for the Hopf bifurcation, given in section [II A[ can be 
obtained by setting 4^ = in the linearized radial equation, 



Since 9 does not appear in the equation for r, the periodic orbits of the system are given by 

r, ,»r 1 V- i t \ t f2Nr\ [2Nr\ 1 

F r {Nr,u),5) = - ^ a m , .{8,w)J m I — \ J m+1 I — I = -y-^ . 



m— 0,oo 



(32) 



These curves are plotted in fig. Ufor cj — 2 and 7 = 0.01, 0.001, 0.0001 and various values of 8. Corresponding to 
these oscillations, the cavity field amplitude oscillates with frequency Q + Fi (Nr,uj,8) and amplitude ey / 2Nr \F\: 

^Leading oscillatory term in \a\ 2 ^ = 2Nre 2 \F (Nr, uj, 8)\ cos ((a) + F z (Nr, lo, 8)) t + Q) , (33) 
where £ is a constant. 

Although the equation governing the periodic orbits is simple, the multiple ranges of the function F r (Nr,u),$), 
whose contours are plotted in fig. [6] for ui = 2, result in multi stability. Its turning points define the positions of the 
saddle node bifurcations, which map out the number of periodic orbits existing in parameter space, as shown in fig. 
[2j (The curves shown in fig. [2] were calculated using 10 terms in the sum.) Expanded in a Taylor series as a function 
of r 2 about zero; 

F r (Nr, R, u) = F r0 (ft, u) + r 2 F rl (ft, u) + r 4 F r2 {ft, u) + ■ ■ ■ . (34) 
the linear term; F r \ (ft,u)), defines the criticality of the Hopf bifurcation. The Hopf bifurcation is super-critical, 
creating a stable periodic orbit, if F r i (ft,u) < 0, which is the case here for 7 small if 8 > — yj 8 "^+ l . 

For larger values of oj the oscillations occur at radii with greater values of N \A\ (= Nr). Fig. [S] compares the 
amplitudes |oj and N \ A\ with the bifurcation diagram for u) = 10 and 7 = 0.00001. For instance in the experiment 
described in [22| the upper bound for the magnitude of epsilon implies that oscillatory behaviour occurs for N on the 
order of 10 and multi-stable behaviour occurs for N on the order of 500. 

Here we will not consider the case with 8 > 0, which corresponds to red detuning, except to note that the dynamics 
is more complicated and deserves a separate study. While periodic orbits, similar to those discussed here, exist there 
are other orbits as well, associated with the Hopf bifurcation, and many of these undergo period doubling (see fig. [2]) 
to chaos. 



C. N nonidentical mechanical resonators and synchronization 



If the frequencies and or damping of each individual mechanical resonator differs, reduction to a single collective 
variables is no longer possible. However the results of the previous section can be generalized to give a set of N 
coupled amplitude equations. Here we consider the case where the linear frequency of the mechanical resonators are 
approximately the same; w, = to + "fAuji. The equations of motion, (I14p . then become 

= - (1 + i5) a - ia V N t X t - it , 
at ^— ' . , 

Xi = - (uj 2 + 2wAw i ) X t - -L-l \ a \ 2 - 2 1 X l . 
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7(NG) e >/(NG) e 



FIG. 4: The amplitudes (N \A\ = Nr) of the periodic orbits of the system calculated from the amplitude equations as a function 
of VNG e for cj = 2 and 7 = 0.01, 0.001, 0.0001 and various values of 5. In a) 5 = -2, b) 5 = -5., c) 5 = -9, d) 6 = -10. The 
unstable periodic orbits are given by dashed lines and the stable one are given by solid lines. 




FIG. 5: The amplitudes \a\ and N \A\ as a e is increased for u = 10 and 7 = 0.00001 compared with the saddle node bifurcations 
that create them. Figure part a) is the cavity amplitude \a\, b) iV|v4| and c) the bifurcation diagram. 

As in the previous section amplitude equations, as function of a slow time (r), can be derived for the dominant 
oscillatory term [ll[ 

X l = X Q + (A t (T)e iult + c.c.) = X Q + 2 \A, (t)| cos(wt + 6 t ) . (36) 
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FIG. 6: The contours of the function F r (Nr,u,S) for uj — 2 plotted as a function of (Nr,S), calculated using 10 terms in the 
sum of products of Bessel functions. Stable periodic orbits exist in green shaded regions. There are no periodic orbits in the 
regions enclosed by red lines where F r (N \A\ , 2, 5) < 0. 



Taking a sum as before and rewriting this as one oscillatory term, 



X = i NiXi = X Q + i Mr)e iu;t + c.c.^J = X + 2 \A\ cos (wi 



(37) 



we can see that \A\ = r acts as a dynamical order parameter for the mechanical resonators, in the sense of Kuramoto 



(38) 



As before, we can use Bessel functions to work out the cavity amplitude response and substitute this back into the 
equations for the individual oscillators to give the amplitude equations, 



dr ' v ' V ^ 

m— 0,oo x 

where a m (S, uA = a mr (S, uA + \a m i (5, ui\ is defined in the previous section. 



UJ 



(39) 
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1. Two sets of nonidentical mechanical resonators 



In terms of two sets of oscillators this becomes 



dr 



(1 + iAwi) A 1 + GNil 2 (A 1 + A 2 ) F (\Ai + A 2 \) , 

(40) 



^ = - (1 + iAw 2 ) A 2 + GN 2 e 2 {Ax + A 2 ) F (\Ax + A 2 \) . 
dr 

If the Aoji are equal they do not effect the radial motion and we still have 

= -r + Ge 2 NrF r (Nr, w, S) , (41) 

dr 

which implies that N 2 r 2 = r 2 +r 2 + 2r\r 2 cos(6 2 — 9\) is a constant of the motion. Substituting this into the equations 
for Ai results in a linear system whose symmetrical solution N\A 2 = N 2 Ai is stable. So apart from some transients 
the individual oscillators synchronize, d ( NlA2 ~ N2Al *> = _ ^ _). jAw) (N\A 2 — N 2 Ai), as noted before. 

If the Awi are not equal the dynamics of the system, which is a function of the relative phase 4> = 9 2 — 9\ only, is 
given by the nonlinear system 

^P- = -n + l 2 GN x (nF r (Nr) + r 2 (F r (Nr) cos 4>- F, (Nr) sin </>)) , 
dr 

^ = -r 2 + l 2 GN 2 (r 2 F r (Nr) + n (F r (Nr) cos <f> + F t (Nr) sin 0)) , (42) 

^ A - . -2„jp /AT UfAT AM , ( N ^l N l T A \ , ,,r^f N 2n , N^ 2 



dr = Aw 2i + e z GF t (Nr) UN 2 - Ni) + I - j cos cj>\ + l A GF r (Nr) \ ~T^ + ) sin< 



where F itT (Nr) = F iiT (Nr,u>,S), Nr = \Ai + A 2 \ = yjr\ + r\ + 2r x r 2 cos cj) and Acj 2 i = Aw 2 - Acji. For Ni = N 2 
we can assume that Auj 2 i > as the transformation (AW21 — > — Auj 2 \, <j) — > —<fi) and (r*i — > r 2 and vice versa) leaves 
the equations unchanged. The coupling, however is strong rather than weak and the system cannot be reduced to a 
phase model. But it is nevertheless useful to compare our results with those of similar phase and phase amplitude 
models, @, |, HO, H, H3 . 

In the simplest two-oscillator phase model (4> = Alu — Ks'mtfi with = 9 2 — 9%) there are two critical points, 
approximately an ina-phase and an out-of-phase solution. One of the critical points is stable, for |Aw| sufficiently 
small (| Aw | < K). Unsynchronized motion occurs when the critical points are lost via a saddle node bifurcation 
(| Aw | > K). More complex models include a sin 20 term in which case the in-phase solution (<p ~ 0) may loose 
stability to a stable out-of-phase solution (<f) 7t). The model here can also be discussed in terms of the stabilities of 
in-phase and out-of-phase solutions. However the 'unsynchronized behaviour' occurs as a transient state, resembling 
the transient rotational motion of a damped nonlinear pendulum started near to the separatrix of the undamped 
system. Similar motion has been noted for other systems with multi stability [23 |. 

Nonzero AW21 breaks the symmetry and the in-phase critical points, which are still stable for AW21 very small, 
only exist with r\ ^ r 2 . Their relative sizes as |Aw2i| is varied in shown in fig. [7]b). As |Aw2i| is increased they 
loose stability via a Hopf bifurcation, fig. [7] a). This creates a stable periodic orbit which does not initially envelope 
the origin. However, in a bifurcation scenario typical of large amplitude coupling 2], it grows rapidly to enclose the 
origin. (In (ri,r 2 ,<fi) space this transition is a heteroclinic bifurcation with saddles at r\ m2 = 0, 4> = ±f ■) Transient 
unsynchronized motion results for solutions started near the (unstable) in-phase solution, where solutions appear 
unbounded in-phase, but eventually become trapped by a stable out-of-phase solution. (In fact the out-of-phase 
solutions are only unstable for AW21 very small, where they exist at large amplitude, fig. He)). The bifurcation 

diagram fig. [7] was created using the package MatCont with F (\A\ 2 , k, tSj approximated by the first four terms in 

its Maclaurin series in \A\ 2 , Ni = N 2 , 7 = 0.0001, w = 2 and «5 = -1.5. 

If we consider only the solutions started near the in-phase solution then, for sufficiently large y/N\Ge > 0.5, 
increasing |Aw2i| engenders a loss of synchronization, see fig. [8] A heteroclinic bifurcation provides the real boundary 
for loss of synchronization and eventually solutions synchronize into an out-of-phase solution. In the unsynchronized 
behaviour the radii execute fairly large oscillations. However the oscillations in the cavity amplitude are not large. A 
typical example is shown in fig. |S]for w = 2, 6 = —1.5, 7 = 0.0001, \fN\Ge = 2 and AW21 = 0.04 starting near the 
unstable in-phase solution. 

If the N are not equal the bifurcation diagram is not symmetrical in AW21. However apart from this it is not 
dissimilar. The in-phase solution with = occurs for £1 = 2£l and looses stability as |Aw2i| is increased away from 
zero, eventually stabilizing on an out of phase solution. 
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FIG. 7: The bifurcation diagram for 2 mechanical resonators for ATi = N2, 7 = 0.0001, uj — 2 and 5 = —1.5. The in-phase 
solutions are stable outside the shaded cyan regions. The out-of-phase solutions are stable outside the slice near the horizontal 
axis. They are singular at Aoj2i = and unstable for |Aw2i| small, where they occur for very large values of n. In the 
unshaded regions both in-phase and out-of-phase solutions are stable, but have different basins of attraction, b) shows the 
in-phase solution in (n, r2) space as ACJ21 is varied, n +T2 remains approximately constant, c) shows the out of phase solution 
in (n, r-2) space as Aaj2i is varied. 
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FIG. 8: Transient unsynchronized motion for 2 nonidentical mechanical resonators for Ni = N2, oj = 2, 5 = —1.5, 7 = 
0.0001, y/N\Ge = 2 and Aoj2i = 0.04. Started near the in-phase solution, in the blue shaded region of fig. [7] the transient 
unsynchronized motion is only temporary. Eventually solutions are trapped by an out-of-phase solution ((f) mod (27r) — > 7r). 
The variables are plotted against time. In a) n and V2 are shown in solid lines, the collective variable r is dashed and cjf>mod(27r) 
is dotted, b) is a plot of <j> and c) is a plot of the amplitude of GN \ a\ 2 . 



2. Three sets of nonidentical mechanical resonators 
The system for N sets of mechanical resonators 



where (Nr) 



dr " 



-n + € 2 n,g 



Au, + e 2 N t G 



r s (F r {Nr) cos (6j - 9$ - F { (Nr) sin [dj -9,)) 

3 

Y, — {Fi (Nr) cos (9 3 - 9,) + F r (Nr) sin (dj - 9,)) 



(43) 



9j), may be reduced to 2N — 1 equations of motion because the 



equations above are only functions of the relative phase: 4n = —6{. So three mechanical resonators are described 
by five equations of motion for n , r%, rs, (f>i, (f>2- If the Aw^ are equal the model can be reduced to that for a single 
collective variable. In fact if any two of the Awi are equal then those two resonators can be thought of as one. Using 
the notation Aw^- = Aui — Auij, the three oscillator case reduces to the two oscillator case if A1J21 — or if AW32 = 
or if AW21 + AU32 = 0. Fig. O shows a typical example of loss of synchronization for AU21 + AW32 small. 

From a dynamical point of view the three resonator case has only one in-phase motion (fa w 0) and one out-of- 
phase motion with (f>i w tt, fa ~ or the other way round. (The case with both fa n is dynamically the same 
as fa w 7T, fa 0.) So as before we can think in terms of the in-phase and out-of-phase solutions. (This is not the 
case for N > 4.) Otherwise the bifurcation diagram is more complicated involving two sets of Hopf curves, however 
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FIG. 9: Transient unsynchronized motion of the in-phase solution for 3 sets of mechanical nonidentical resonators for Ni equal, 
uj — 2, 8 — —1.5, 7 = 0.0001, y/NiGe = 2, Aoj2i = 0.04 and Au2i = 0.045. Eventually solutions are trapped by an out-of-phase 
solution (<f>i mod (2-k) — > n here). The variables are plotted against time. In a) rj are shown in solid lines, the collective variable 
r is dashed and </>imod (2tt) are dotted, b) is a plot of 4>i and c) is a plot of the amplitude of GN |a| 2 . 



if AW21 and 5Ail>32 are close the Hopf curves are also close. In contrast if they differ, as shown in fig. [TU] where we 
have taken AW21 = 5Aw32, three unstable regions result. The most complex behaviour occurs in the blue region in 
which the in-phase solution is unstable to both <f>i and the motion may switch from librational to rotational motion 
in one or both of the <f>i apparently randomly. 



III. QUANTUM MECHANICAL DESCRIPTION 
A. Quantum mechanical model 

The classical model derived in section U and analysed in section [TT| is given by the classical Hamiltonian 



N 



N 



H = H6 \a\ 2 + hoji \a\ 2 + h (e*a + ea*) + figi \a\ 2 Xi 



(44) 



where the microwave cavity amplitude a — x c + iy Cl and the dimensionless canonical positions x c and Xi, and their 
conjugate momenta y c and yi respectively, satisfy the Poisson bracket relations 



{xc,Vc} - 2h , 
{xi, yj} — 2^ij ■ 



(45) 
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0.5 1 1.5 2 Z5 3 3.5 A 
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FIG. 10: The bifurcation diagram for three mechanical resonators for Ni equal, uj — 2, S — —1.5, 7 = 0.0001, AtJ2i = 5Aoj32. 
The in-phase solution is stable outside the shaded cyan and blue regions. The out-of-phase solutions are stable outside the grey 
regions. In the blue region the in-phase solution is unstable to both cj>\ and <f>2 and the motion is transiently unsynchronized, 
eventually setting on an out-of-phase solution. Once again the out-of-phase solutions tend to exist where at least some of the 
n take larger values. 



y/2E c L x c and qi 



The original canonical positions $ 
and pi — \/2mE m y$ respectively, satisfy the canonical Poisson bracket commutation relations 



2E,„ 
muj- 



and their conjugate momenta Q — \/2E c Cq y c 



{<f>,Q} =2h{x c 



Vc) = 1 , 



(46) 



The quantum mechanical description of the Hamiltonian dynamics matches that obtained by canonical quantisation 
of the classical Hamiltonian. We promote the canonical position and momenta <P and Q of the microwave cavity to 
the quantum mechanical operators 4> and Q respectively. Similarly, we promote the nano-mechanical resonator 
positions qi and momenta pi to the quantum mechanical operators and respectively. We then define the 
annihilation operators for the microwave cavity field mode a, and the nano-mechanical vibrational modes hi. We 



again have the dimensionless microwave cavity quadrature operators x c — 



the dimensionless nano-mechanical positions and momenta cti 
We have the commutation relations for the quantum operators 



and y i 



and y c = — i| 

-ii (i 



and 



it 



h.bj 



= 1, 



= Siji 



[x c ,y c 



ihSijl . 



respectively. 



(47) 
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in terms of which the corresponding quantum Hamiltonian in the interaction picture is given by 



U = HS I a f a + ^ 



1 



JV 



2J 



bik + ~ 



N 



h ( e*a + ea^ ) + 



(48) 



For a realistic device we adopt a dissipative model. We model both the microwave cavity resonator and the mechanical 
resonators as being damped in zero temperature heat baths. A zero temperature heat bath for the cavity is certainly 
justified as the typical microwave cavity is at mK temperature and thus very close to zero [28| . The N zero temperature 
heat baths for the N nano-mechanical resonators are not as good an approximation. However, the mean thermal 
occupation of the zth bath rij ^ does not enter the semi-classical equations, and thus the semi-classical bifurcation 
structure studied in section |TT] is the correct one. The amplitude decay for the microwave cavity is n, and for the 
ith nano-mechanical resonator is 7.;. We then describe the dissipative dynamics with the master equation (with weak 
damping and the rotating wave approximation for the system-environment couplings) 



dp 
dt 



H,p 



N 



t t 



a) ap — pa)a \ + 7, febipbi 

i=l 



bi hip - ph bA , 



(49) 



where p is the density matrix of the coupled system. 

Corresponding to the the classical description, we are interested in the M collective quantities Xi and Y i defined 

by 



Y, 



-E 



gjXj 



9jVj 



(50) 



We can define creation and annihilation operators for these collective mechanical modes, 

t 1 V ;t 



(51) 



Bi 



jeSi 



where from (147[) . we can show that the commutation relations for the new collective operators are 



Gi 



5i.il . 



(52) 



B. Fokker-Planck like equation 

From the Master equation (|49p , we proceed by deriving a Fokker-Planck like equation for the nano-electromechanical 
system which is the equation of motion of the Positive P function P(x)- The Positive P function is the Fourier 
Transform of the expectation of the normally-ordered characteristic function, 

P W = ^ 2M +2 J (e iX2M+2§M \ iX2M+lilM ■ ■ ■ e iA ^i t e iA3jll e iA2at e iAl<i ^ e^ A '*dA , (53) 

where 

X = \ a (3 pi vi ^2 v 2 ■ ■ ■ pm vm] T , , , 

T (54) 

A = [ Ai A2 • • • X2M+2 \ 

We follow the procedure outlined in [19| . Using the appropriate commutation relations, we arrive at the Fokker-Planck 
like equation 

^ = "E 4 [AMI, PW + \ E [B«n s P(X) , (55, 
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where the drift term vector A(x) is 



A(x) 



-ie - (« + iS) a - igaj^j^i (Mi + v i) 
ie - (k - W) /3 + i±/3 E<=i ^ (Mi + 

- (7 + iwi) fii - i^-a/3 

— (7 — ilui) v\ + 



(7 + ioj M ) MM -i%-a/3 
(7 - iwAf) km + 



(56) 



and the diffusion term matrix B(x)B(x) T is 



B(X)B( X ) 5 









-if a 


• • 


• -i^fa 














if/3 •• 





if*-/3 


-if Q 








■ • 











if/3 





•• 
















• • 











if*-/3 





•■ 









(57) 



If we consider only the drift term of the Fokker- Planck like equation (f5"5j) . and make the mappings j3 h-> a* and 
u 1 y v* to reduce the phase space dimension by half onto the semi-classical phase space (the positive P function has 
twice the dimensionality of the classical phase space), then we obtain the semi-classical equations of motion. 



C. Quantum spectra 



A future direction for research, that builds on the work of this paper, is the investigation of the quantum physics 
associated with the multi-stable semi-classical limit cycles. As a starting point, in this section we calculate the 
linearised spectrum as we increase the driving strength to approach the first Hopf bifurcation at the supercritical 
Hopf line for blue detuning (S < 0) in Fig. [2] and Fig. [3] We do this calculation for the case of a single group 
of nano-mechanical resonators following the procedure of [29j . For a single group, using the dimensionless notation 
where we have re-scaled the coupling coefficients and time by the cavity dissipation rate k, we have the stochastic 
differential equations of motion corresponding to the Fokker-Planck like equation (|5"5"|) 



^=A(*)+B( X )E(f), 



where 



X = [ a /3 n v ] T , 



(58) 



(59) 



the drift term vector A(x) is 

Mx) = 

the diffusion term matrix B(x)B(x) T is 

B(x)B( X ) T 



(1 + iS) a - i^aN (h + l>)- ie 
(1 - iS) /3 + iJjSJV (p + v)+ie 

- (7 + iu)fj, - ij-a(3 

— (7 — icj) v + ifa/3 



-ifa 

if/3 

-ifa 

if/3 



(60) 



(61) 
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and E(t) is the noise process. The principal matrix square root of the diffusion matrix B(x)B(x) J 



is 



B( X )=B( X )* = ^ 



y/a — \y[a 

yf$ iy/P 

-iy/a yfa 

iy/P 







(62) 



The diffusion matrix and its square root have determinants 



det{B(x)B( X ) T } = -GW 
det {B( X )} = \G 2 aP, 



(63) 



and the two matrices are thus positive definite on the semi-classical manifold where (3 = a*. We see that the off- 
diagonal terms with the factors of i in the matrix square root B(x) will take the solution off the semi-classical manifold, 
and will lead to quantum correlations. 

We will linearise these equations of motion about the semi-classical fixed points we obtained in section [TTJ In terms 
of the stochastic differential equations above, we make the mappings /3 i— > a* and u i— > v* to reduce the phase space 
dimension by half onto the semi-classical phase space (the positive P function has twice the dimensionality of the 
classical phase space). The linearised stochastic differential equations are then 



where our Jacobian matrix M is 
df(Xo) 



M 



- (1 + i5) - i§JV (a»o + A*o) 



"if "5 



D'E(t) 







(64) 



(l-iS)+iiiVOuo + Mo) 
-iyao 

i T a n 



-ihao 


-ioQ!0 


■ i * 


• 1 * 


(7 + id) 










- (7 - iw) 

Xq = 5(^0 + Mo) an d our diffusion matrix about the semi-classical fixed points D = B(x )B(x ) T is 



(65) 



D 






-if 



-ifa 




1^05 



(66) 



ifa£ 

The linearised normally ordered moments at steady state can be expressed in terms of these matrices [l9| 



1 

2^ 



dr = — (ii7I - M)" 1 D (-im 
2vr V ; V 



(67) 



We plot the microwave cavity component of these quantum noise spectra in Fig. [TT] We see that in Fig. [TT] (a) as 
the Hopf bifurcation is approached, the spectrum becomes more sharply peaked at two frequencies. The frequency 
corresponding to the Hopf bifurcation — the magnitude of the two purely imaginary eigenvalues — is the peak at 
the mechanical frequency uj. The second, shorter but broader peak, is at the detuning 8. For a drive detuned exactly 
on a sideband, these two peaks coincide. Beyond the supercritical Hopf bifurcation, the semi-classical fixed point 
is no longer stable and we enter the regime dominated by the first stable limit cycle, where we have the oscillatory 
motion analysed by semi-classical amplitude equations in section |Hj However, we can continue to linearise about this 
point, and show the results in Fig. [TT](b). The two peaks begin to converge as the driving strength and coupling are 
increased. 

The spectra calculated here correspond to the stationary fluctuations in the cavity field. The power spectrum of 
these fluctuations can be directly measured by homodyne detection. Below the Hopf bifurcation, the noise power 
spectrum of these fluctuations is peaked at a frequency associated with the decay of fluctuations back to the fixed 
point. The width of the peaks gives the time scale of this decay. Above the Hopf bifurcation, the fluctuations 
decay onto the limit cycles. In our model, there are no thermal fluctuations and all fluctuations are due to intrinsic 
quantum noise manifest as off-diagonal components in the diffusion matrix in (|66[) . A more careful study is required to 
determine if the multiple peak structure evident in Fig. Illf b) is evidence for dissipative quantum switching between 
limit cycles. 
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FIG. 11: Linearised quantum noise spectrum of the microwave cavity Sn(O): (a) approaching the Hopf bifurcation; and (b) 
continuing the linearisation beyond the Hopf bifurcation. The magnitude of the normally ordered cavity spectrum at steady 
state /T^o e ~ lflT ( a (t) a (t + T ) T ) t ^ e nrs t diagonal element of S(£2), is plotted at the frequency fl for varying driving 

amplitude e. Here we have set oj = 10, 5 — —4, 7 = 0.001, N = 1, G = 1, for which the Hopf bifurcation occurs at a driving 
strength of eh « 1.76. 



IV. DISCUSSION AND CONCLUSION 

We have discussed the situation in which multiple mechanical resonators are coupled to a single mode of the 
electromagnetic field in a microwave superconducting cavity. This interaction results in an all-to-all coupling between 
each of the mechanical resonators that is highly nonlinear. However if the oscillators are identical, they synchronize 
and a collective variable can be used to understand the dynamics. Analysis of the dynamics of this collective variable 
(see the bifurcation diagram in fig. [5]) reveals the prevalence of periodic behaviour and suggests the use of amplitude 
equations (|2"5|) to describe the dominant oscillation. Even though the amplitude equations involve elliptic functions 
their overall form is relatively simple for small mechanical damping and from them we are able to gain considerable 
insight into the dynamics of the collective variable. 

The form of the amplitude equations imply the presence of multiple periodic orbits and hysteresis (at the bifurcations 
of the periodic orbits). Specific results are also easy to extract; for instance we are able to plot the amplitudes of the 
mechanical resonators as a function of the external forcing for specific values of the other parameters (fig. 13} and to 
locate the saddle node bifurcations of periodic orbits where hysteresis would occur as a result of a slight change in 
the mechanical forcing (fig. fig. 02 and fig. [5]). 

The simplicity of the amplitude equations means that it is straightforward to extend the identical mechanical 
resonator case to one with distinct subgroups of identical oscillators. Considering two and three frequency subgroups 
we are able to give bifurcation diagrams showing the regions where synchronization occurs. Synchronization is lost via 
a mechanism involving a Hopf and heteroclinic bifurcation similar to that found in large amplitude forcing, rather than 
the sniper bifurcation that is involved in small amplitude forcing and, although in a reduced form, in Kuramoto's phase 
model. In spite of this difference there is a single collective variable Nr that functions as a measure of synchronized 
behaviour and that is related to a measurable quantity, the cavity amplitude. 

Given the current interest in fabricating nanomechanical resonators in microwave cavities, our model offers a 
realizable and very controllable way to study synchronization in a system with all-to-all coupling via a common field 
mode. While the equations for our model cannot be reduced to a simple phase model, it offers some advantages over 
more complex naturally occurring examples of synchronization. A particularly important feature is that the measured 
quantity — the cavity field leaving the microwave resonator — has an amplitude that is directly proportional to a 
collective parameter similar to the order parameter introduced in previous studies of synchronization. The need to 
use very low temperatures required for superconducting circuits may seem a disadvantage but in fact leads to a huge 
reduction in noise both for the mechanics and the microwave field. This should lead to especially clean observation 
of multi stability and perhaps even controlled switching between limit cycles. In the long run it also motivates us to 
study the effect of quantum noise on synchronization, and to look for quantum signatures of synchronization which 
will be the subject of a future paper. 
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Appendix A: Experimental parameters 



A number of experiments are described by the model examined in this paper, 



dp i 

dr~ "~ ~ft 



N 

k \2apaJ - a)ap - paJaA + ^ 7, (%bipbi - bi hip - pb t b^j , 



(Al) 



where 



/ 1\ N { \\ N 1 / 1\ 

U = HS (Va + - J + ^ fiwj ^Sj&i + - J + ft (e*a + ea^ + ^ \tfa + - J (fc + k) , (A2) 



and 



&i,6j =5ijl. 



(A3) 



A summary of the different values of the parameters for a selection of these experiments is given in table [I] In terms 
of the dimensionless parameters introduced in section [TT| these become those listed in table |TIJ Note that the detuning 
6 is typically set to be on a mechanical frequency sideband, such that 5 = ujf, and thus while not an experimental 
limitation, the range of S we list in the table is S < Wj. Also note that the maximum driving |e| indicates the maximum 
driving before the cavity becomes nonlinear causing our model to fail. Finally, also note that the factors of 2 in front 
of k and 7^ in table U are present because our k and 7, (as defined by the master equation above) are amplitude decay 
rates not occupation number decay rates. 



Experiment 
Ref. Type 
S 
S 

s 



[30] 
[22] 

[14] 



[U] 
[31] 
[23] 
[15] 
[16] 

[16] 
[16] 



S 
S 

s 

M 
T 

T 
T 



7.49 x 10 9 
5.22 x 10 9 

4 x 10 9 
-j- 10 x 10 9 
7.55 x 10 9 
~ 5 x 10 9 
7.64 x 10 9 
282 x 10 12 



Mode a 

2£ [Hz] 
< 2.88 x 10 6 
230 x 10 3 
400 x 10 3 
-> 1 x 10 6 
302 x 10 3 
490 x 10 3 
382 x 10 3 
4.07 x 10 6 
> 4.9 x 10 6 

50 x 10 6 
50 x 10 6 



Mode bi 



M 

2tt 



Hzl 



< 2.145 x 10 9 



1.04 x 10 6 
1.53 x 10 6 
0.1 x 10 6 
->• 6 x 10 6 
1.41 x 10 6 
2.3 x 10 6 

< 2.434 x 10 9 67 x 10 6 
134 x 10 3 

6.5 x 10 6 
-> 16 x 10 6 
10.74 x 10 6 

8 x 10 6 



2^ P Z ] 
0.67 

< 5.08 
< 1 

^< 6 

< 371.1 
19.2 
248.1 
0.122 

65 
-> 1600 
202.64 
200 



Coupling 

ft [Hz] 
866.7 x 10- 
190.7 x lO - 



49.55 x 10" 
25.03 
27.8 



147.3 
55.6 



TABLE I: Raw experimental coupling values for various systems. The 'Type' column indicates the experimental context: 'S' 
indicates a superconducting microwave coplanar waveguide resonator (a) coupled to a nano-mechanical resonator (bi); 'M' 
indicates an optical cavity (a) coupled to a micro-mechanical membrane (bi); "I" indicates a toroidal micro-resonator (a) 
coupled to a nano-mechanical string resonator (bi); and 'C indicates an opto-mechanical crystal array where an optical mode 
of a cell (a) is coupled to a mechanical mode of a cell (bi). 
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Experiment 


Mode a 




Mode hi 


Coupling 


Ref. 


Type 


|5'| = J|1 |e'| = ^l 


/ U)i 
U i = -f 


i ~/j 
7» = « 


9'i=*t 


[30] 


S 


< 0.722 


0.722 


2.33 x 10~ 7 


6.02 x 10" 


[22] 


s 


< 13.26 < 1.87 x 10 4 


13.26 


2.209 x 10" 5 


1.66 x 10" 


M 


s 


^ 0.5 


0.5 


2.5 x 10" 6 






12 


-> 12 


-46x 10" 6 




M 


s 


< 9.34 


9.34 


1.23 x 10" 3 




[31] 


s 


< 9.39 


9.39 


3.92 x 10" 5 


2.02 x 10" 


[23] 


s 


< 350.8 < 1.28 x 10 4 


350.8 


6.5 x 10" 4 


1.31 x 10" 


[15] 


M 


< 0.066 


0.066 


3 x 10" 8 


1.37 x 10" 


[16] 


T 


< 2.65 


< 2.65 


< 1.33 x 10" 5 








->< 6.53 


->< 6.53 


-»< 3.27 x 10" 4 




[16] 


T 


< 0.43 


0.43 


4.05 x 10" 6 


5.89 x 10" 


[16] 


T 


< 0.32 


0.32 


4 x 10" 6 


2.22 x 10" 



TABLE II: Dimensionless experimental coupling values for various systems. The 'Type' column indicates the experimental 
context: 'S' indicates a superconducting microwave coplanar waveguide resonator (a) coupled to a nano-mechanical resonator 
(bi); 'M' indicates an optical cavity (a) coupled to a micro-mechanical membrane (bi); 'T' indicates a toroidal micro-resonator 
(a) coupled to a nano-mechanical string resonator (bi); and 'C indicates an opto-mechanical crystal array where an optical 
mode of a cell (a) is coupled to a mechanical mode of a cell (bi). 
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